home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / swagg_m.zip / MATH.SWG / 0045_Calculate PI.pas < prev    next >
Pascal/Delphi Source File  |  1993-11-02  |  2KB  |  70 lines
  1. {
  2. LOU DUCHEZ
  3.  
  4. ATTENTION, whoever was trying to calculate PI!  Here's a swell program,
  5. as a follow-up to a recent post of mine about approximating techniques!
  6.  
  7. }
  8.  
  9. program calcpi;  { Calculates pi by getting the area of one-quarter of a
  10.                    circle of radius 1, and then multiplying by 4.  The area
  11.                    is an approximation, derived by Simpson's method: see
  12.                    previous post for explanation of that technique. }
  13.  
  14. uses
  15.   crt;
  16.  
  17. const
  18.   lowerbound  = 0;  { The interval we're evaluating is from 0 to 1. }
  19.   higherbound = 1;  { I put the 0 and 1 here for clarity. }
  20.  
  21. var
  22.   incs    : word;
  23.   quartpi,
  24.   h, x    : real;
  25.  
  26. function y(x : real) : real;  { Feed it an x-value, and it tells you the }
  27. begin                         { corresponding y-value on the unit circle. }
  28.   y := sqrt(1 - x * x);       { A no-brainer. }
  29. end;
  30.  
  31. begin
  32.   { I leave you to do the error-checking on input. }
  33.   clrscr;
  34.   write('Enter a WORD (1 - 32767) for the number of parabolas to do: ');
  35.   readln(incs);
  36.  
  37.   { The answer for a quarter of pi will be accumulated into QuartPi. }
  38.   quartpi := 0;
  39.  
  40.   { H is the interval to increment on.  X is the "middle" x value for each
  41.     parabola in Simpson's method.  Here it is set equal to one interval
  42.     above the lower bound: Simpson's method looks at points on either side
  43.     of "X", so my reasoning is obvious.  Note also that, by magical
  44.     coincidence, the last evaluation will have "X" equal to the higher
  45.     bound of the interval minus H. }
  46.  
  47.   h := (higherbound - lowerbound) / (1 + 2 * incs);
  48.   x := lowerbound + h;
  49.  
  50.   { This loop accumulates a value for pi/4. }
  51.   while incs > 0 do
  52.   begin
  53.     if x < 0 then
  54.       x := 0;
  55.     quartpi := quartpi + y(x - h) + 4 * y(x) + y(x + h);
  56.  
  57.     { Move X two increments to the right, and decrease the number of parabolas
  58.       we still have to do. }
  59.     x := x + 2 * h;
  60.     dec(incs);
  61.   end;
  62.  
  63.   { Simpson's method has you multiply the sum by H/3. }
  64.   quartpi := h * quartpi / 3;
  65.  
  66.   { Print answer. }
  67.   writeln(4 * quartpi : 12 : 8);
  68.   writeln('This has been a display of Simpson''s method.  D''ohh!');
  69. end.
  70.